Raw data
- Study aim: Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data
- Study summary: Sequencing Quality Control Consortium generated two datasets from two reference RNA samples in order to evaluate transcriptome profiling by next-generation sequencing technology. Each sample contains one of the reference RNA source and a set of synthetic RNAs from the External RNA Control Consortium (ERCC) at known concentrations. Group A contains 5 replicates of the Strategene Universal Human Reference RNA (UHRR), which is composed of total RNA from 10 human cell lines, with 2% by volume of ERCC mix 1. Group B includes 5 replicate samples of the Ambion Human Brain Reference RNA (HBRR) with 2% by volume of ERCC mix 2. The ERCC spike-in control is a mixture of 92 synthetic polyadenylated oligonucleotides of 250-2000 nucleotides long that are meant to resemble human transcripts.
- Sample: A(replicates of the Strategene Universal Human Reference RNA (UHRR)), B(replicate samples of the Ambion Human Brain Reference RNA (HBRR) with 2% by volume of ERCC mix 2)
Loading packages
library(data.table)
library(rmarkdown)
library(AnnotationHub)
library(tidyverse)
library(tximport)
library(ggplot2)
library(DESeq2)
library(pheatmap)
library(gridExtra)
library(ggplotify)
library(ggrepel)
library(UpSetR)Setting AnnotationHub
Assign your species of interest
AnnotationSpecies <- "Homo sapiens" # Assign your species
ah <- AnnotationHub(hub=getAnnotationHubOption("URL")) # Bring annotation DBRunning AnnotationHub
ahQuery <- query(ah, c("OrgDb", AnnotationSpecies)) # Filter annotation of interest
if (length(ahQuery) == 1) {
DBName <- names(ahQuery)
} else if (length(ahQuery) > 1) {
DBName <- names(ahQuery)[1]
} else {
print("You don't have a valid DB")
rmarkdown::render()
}
AnnoDb <- ah[[DBName]] # Store into an OrgDb object
# Explore your OrgDb object with following accessors:
# columns(AnnpDb)
# keytypes(AnnoDb)
# keys(AnnoDb, keytype=..)
# select(AnnoDb, keys=.., columns=.., keytype=...)
AnnoKey <- keys(AnnoDb, keytype="ENSEMBLTRANS")
# Note: Annotation has to be done with not genome but transcripts
AnnoDb <- select(AnnoDb,
AnnoKey,
keytype="ENSEMBLTRANS",
columns="SYMBOL")Checking out the AnnotationHub output
# Check if your AnnoDb has been extracted and saved correctely
class(AnnoDb)## [1] "data.frame"head(AnnoDb)## ENSEMBLTRANS SYMBOL
## 1 ENST00000543404 A2MP1
## 2 ENST00000566278 A2MP1
## 3 ENST00000545343 A2MP1
## 4 ENST00000544183 A2MP1
## 5 ENST00000286479 NAT2
## 6 ENST00000520116 NAT2Defining file name and path for .sf files
.sf files have been created from fastq data by salmon
# This code chunk needs to be written by yourself
# Define sample names
SampleNames <- c(paste0("A", 1:3), paste0("B", 1:3))
# Define group level
GroupLevel <- c("A", "B")
# Define contrast for DE analysis
Contrast <- c("Group", "B", "A")
# Define a vector for comparing TPM vs Counts effect
TvC <- c("TPM", "Counts")
# Define .sf file path
sf <- c(paste0("salmon_output/",
SampleNames,
".salmon_quant/quant.sf"))
# Define sample groups
group <- c(rep("A", 3), rep("B", 3))
# Create metadata
metadata <- data.frame(Sample=factor(SampleNames, levels=SampleNames),
Group=factor(group, levels=GroupLevel),
Path=sf)
rownames(metadata) <- SampleNames
# Explore the metadata
print(metadata)## Sample Group Path
## A1 A1 A salmon_output/A1.salmon_quant/quant.sf
## A2 A2 A salmon_output/A2.salmon_quant/quant.sf
## A3 A3 A salmon_output/A3.salmon_quant/quant.sf
## B1 B1 B salmon_output/B1.salmon_quant/quant.sf
## B2 B2 B salmon_output/B2.salmon_quant/quant.sf
## B3 B3 B salmon_output/B3.salmon_quant/quant.sfConverting .sf files to txi list
- txi_tpm: stores TPM with the argument “countsFromAbundance=”lengthScaledTPM"
- txi_counts: stores original counts
- Despite that the TPM matrix is not recommended as an input for DE analysis due to the fact that it doesn’t account gene length, TPM stored in a tximport (by tximport(…, countsFromAbundance=“lengthScaledTPM”)) can be used for DE analysis by being inputted with the DESeqDataSetFromTximport() funtion in DESeq2 workflow as the gene length is automatically adjusted by DESeqDataSetFromTximport().
- In this project, two txi objects were created with or without the countsFromAbundance=“lengthScaledTPM” argument and compared in downstream DE analysis.
- If you don’t want gene-level summarization, set txOut=TRUE.
- References: tximport doc, DESeq2 doc “Why unnormalized counts?”, Soneson et al. 2016, Developer Dr. Love’s comment, Harvard Chan Bioinformatics Core workshop
# Assign sample names to the input (.sf) file path
names(sf) <- SampleNames
# Run tximport
# tpm vs original counts
# input sf: a factor of all .sf files' path
txi_tpm <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
countsFromAbundance="lengthScaledTPM", # Extracts TPM
ignoreTxVersion=T)
txi_counts <- tximport(sf,
type="salmon",
tx2gene=AnnoDb,
ignoreTxVersion=T) Exploring the txi outputs
# tpm
head(txi_tpm$counts)## A1 A2 A3 B1 B2 B3
## A2MP1 11.828289 9.640769 4.3214387 10.768983 9.930742 23.582072
## A3GALT2 0.000000 0.000000 0.9922462 4.934298 8.216794 9.121198
## A4GNT 3.878022 1.983684 0.9962060 5.005617 7.978542 4.061164
## AACSP1 99.547876 123.181812 54.2225647 15.620886 26.531886 22.692137
## AADACL2 0.000000 0.000000 0.0000000 1.016813 0.000000 0.000000
## AADACL4 1.951536 0.000000 0.9857014 0.000000 0.000000 0.000000dim(txi_tpm$counts)## [1] 12202 6# counts
head(txi_counts$counts)## A1 A2 A3 B1 B2 B3
## A2MP1 12.118 14.265 3 10.045 12.188 16.12
## A3GALT2 0.000 0.000 1 5.000 8.000 9.00
## A4GNT 4.000 2.000 1 5.000 8.000 4.00
## AACSP1 97.000 129.000 59 15.000 24.000 24.00
## AADACL2 0.000 0.000 0 1.000 0.000 0.00
## AADACL4 2.000 0.000 1 0.000 0.000 0.00dim(txi_counts$counts)## [1] 12202 6Creating DESeq objects from txi and VST
- Note: The tximport-to-DESeq2 approach uses estimated gene counts from the transcript abundance quantifiers, but not normalized counts.
- The DESeqDataSetFromTximport() function generated an DESeq object (aka dds) with the txi input.
- vst() was run to perform variance stabilizing transformation instead of rlog() which takes longer time with similar characteristics.
- The vsd object created by vst() is used for not DE analysis but QC.
- References: DESeq2 doc “Transcript abundance files”, DESeq2 doc “Variance stabilizing transformation”
# Set a function creating dds and vsd
dds_vsd_fn <- function(txi) {
# Create a DESeq object (so-calledd dds)
des <- DESeqDataSetFromTximport(txi,
colData=metadata,
design=~Group)
# Create a vsd object (so-called vsd)
ves <- vst(des, blind=T)
# Output them as a list
return(list(dds=des, vsd=ves))
}
TPM <- dds_vsd_fn(txi_tpm)
Counts <- dds_vsd_fn(txi_counts)
# Outputs
# dds: TPM/Counts[[1]] or TPM/Counts[['dds']]
# vsd: TPM/Counts[[2]] or TPM/Counts[['vsd']]Exploring created dds
# TPM
TPM[[1]]## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(1): counts
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group Pathhead(counts(TPM[[1]]))## A1 A2 A3 B1 B2 B3
## A2MP1 12 10 4 11 10 24
## A3GALT2 0 0 1 5 8 9
## A4GNT 4 2 1 5 8 4
## AACSP1 100 123 54 16 27 23
## AADACL2 0 0 0 1 0 0
## AADACL4 2 0 1 0 0 0# Counts
Counts[[1]]## class: DESeqDataSet
## dim: 12202 6
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(0):
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group Pathhead(counts(Counts[[1]]))## A1 A2 A3 B1 B2 B3
## A2MP1 12 14 3 10 12 16
## A3GALT2 0 0 1 5 8 9
## A4GNT 4 2 1 5 8 4
## AACSP1 97 129 59 15 24 24
## AADACL2 0 0 0 1 0 0
## AADACL4 2 0 1 0 0 0Exploring created vsd
# TPM
TPM[[2]]## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): A1 A2 ... B2 B3
## colData names(4): Sample Group Path sizeFactor# Counts
Counts[[2]]## class: DESeqTransform
## dim: 12202 6
## metadata(1): version
## assays(1): ''
## rownames(12202): A2MP1 A3GALT2 ... ZXDB ZYXP1
## rowData names(4): baseMean baseVar allZero dispFit
## colnames(6): A1 A2 ... B2 B3
## colData names(3): Sample Group PathEstimating size factors, dispersions, and conducting Wald Test
- Dispersion is calculated as a measure of variation instead of variance since variance gets larger when gene expression gets higher.
- Wald test is the default setting of DESeq2 which tests null hypothesis between two groups. You should use Likelihood ratio test (LRT) when comparing more than two groups.
- Messages when “Counts <- DESeqPrep_fn(Counts)” was run:
using ‘avgTxLength’ from assays(dds), correcting for library size gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates
- References: Harvard Chan Bioinformatics Core workshop I, Harvard Chan Bioinformatics Core workshop II, Harvard Chan Bioinformatics Core workshop III, DESeq2 "Wald test indivisual steps, DESeq2 doc “Likelihood ratio test”
# Set a function estimating size factors, dispersions, and perform wald test
DESeqPrep_fn <- function(List) {
List[[1]] <- estimateSizeFactors(List[[1]])
List[[1]] <- estimateDispersions(List[[1]])
List[[1]] <- nbinomWaldTest(List[[1]])
return(List)
}
# Update dds with the function
Counts <- DESeqPrep_fn(Counts)
TPM <- DESeqPrep_fn(TPM)Exploring size factors
sizeFactors(Counts[[1]])## NULLsizeFactors(TPM[[1]])## A1 A2 A3 B1 B2 B3
## 1.2059202 1.0600776 0.7425445 1.0155827 1.0482637 1.0759300# Size factors don't exist in the Counts dds!
# Normalization factors are calculated in the Counts dds instead!
assays(Counts[[1]])## List of length 6
## names(6): counts avgTxLength normalizationFactors mu H cooksassays(TPM[[1]])## List of length 4
## names(4): counts mu H cookscolData(Counts[[1]])## DataFrame with 6 rows and 3 columns
## Sample Group Path
## <factor> <factor> <character>
## A1 A1 A salmon_output/A1.sal..
## A2 A2 A salmon_output/A2.sal..
## A3 A3 A salmon_output/A3.sal..
## B1 B1 B salmon_output/B1.sal..
## B2 B2 B salmon_output/B2.sal..
## B3 B3 B salmon_output/B3.sal..colData(TPM[[1]])## DataFrame with 6 rows and 4 columns
## Sample Group Path sizeFactor
## <factor> <factor> <character> <numeric>
## A1 A1 A salmon_output/A1.sal.. 1.205920
## A2 A2 A salmon_output/A2.sal.. 1.060078
## A3 A3 A salmon_output/A3.sal.. 0.742545
## B1 B1 B salmon_output/B1.sal.. 1.015583
## B2 B2 B salmon_output/B2.sal.. 1.048264
## B3 B3 B salmon_output/B3.sal.. 1.075930Plotting the size factors in TPM
- The size factors are only available from TPM dds
# Extract and save the size factors in a data frame
sizeFactor <- as.data.frame(round(sizeFactors(TPM[[1]]), 3))
colnames(sizeFactor) <- 'Size_Factor'
sizeFactor <- sizeFactor %>%
rownames_to_column(var="Sample") %>%
inner_join(metadata[, 1:ncol(metadata)-1], by="Sample")
# Create a plot comparing the size factors by sample
ggplot(sizeFactor, aes(x=Sample,
y=Size_Factor,
fill=Group,
label=Size_Factor)) +
geom_bar(stat="identity", width=0.8) +
theme_bw() +
ggtitle("Size Factors in TPM-DESeq") +
geom_text(vjust=1.5) +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) + ylab("Size Factor")
Plotting nornalization factors in the Counts
- DESeq2 performs an internal normalization where geometric mean is calculated for each gene across all samples. The counts for a gene in each sample is then divided by this mean. The median of these ratios in a sample is the size factor for that sample.
- Blue dashed line: normalization factor = 1
- Colored dots: normlization factors per gene (y-axis) in each sample
- Box plots: distribution of the normalization facters in each group (see the second plot)
- Reference: DESeq2 doc “Sample-/gene-dependent normalization factors”
# Extract and normalization factors in a data frame
normf <- as.data.frame(normalizationFactors(Counts[[1]])) %>%
gather(Sample, Normalization_Factor) %>%
inner_join(metadata[, 1:2], by="Sample")
normf$Sample <- factor(normf$Sample, levels=SampleNames)
normf$Group <- factor(normf$Group, levels=GroupLevel)
# Create a scatter plot showing distribution of normalization factors
normFactors_plot <- ggplot(normf,
aes(x=Sample, y=Normalization_Factor)) +
geom_jitter(alpha=0.5, aes(color=Group)) +
# Add a boxplot to provide statistics in each sample
geom_boxplot(aes(x=Sample, y=Normalization_Factor),
outlier.shape=NA, alpha=0.5) +
theme_bw() +
ggtitle("Normalization Factors in Counts-DESeq") +
theme(axis.text.x=element_text(angle=45,
vjust=0.5)) +
ylab("Normalization Factor / Gene") +
# Add a dashed horizontal line to indicate where normalization factor=1
geom_hline(yintercept=1,
color="blue",
linetype="dashed")
# Print the normalization factor scatter plot
print(normFactors_plot)
# Print the same plot with larger y-magnification in order to observe the box plot
normFactors_plot +
ylim(0.5, 1.5)
Results & Discussion 1:
- Patterns of the size factors (from TPM) and the normalization factors per gene (from Counts) by sample were similar in terms of magnitude
Setting functions for QC plots
- Reference: DESeq2 doc “Principal component plot of the samples”, DESeq2 doc “Heatmap of the sample-to-sample distances”
# Assigne what to compare
GroupOfInterest <- Contrast[1]
# Set a function for a PCA plot
QCPCA_fn <- function(inputList, Title) {
plotPCA(inputList[[2]], # takes vsd
intgroup=GroupOfInterest) + theme_bw() + ggtitle(Title)
}
# Set heatmap annotation
ColOfInterest <- !colnames(metadata) %in% c("Sample", "Path")
HeatmapAnno <- as.data.frame(metadata[, ColOfInterest])
rownames(HeatmapAnno) <- SampleNames
colnames(HeatmapAnno) <- colnames(metadata)[ColOfInterest]
# Set a function for a correlation heatmap
QCcorrHeatmap_fn <- function(inputList, Title) {
# Extract transformed count matrix
mtx <- assay(inputList[[2]]) # takes vsd
# Calculate correlation and store in the matrix
mtx <- cor(mtx)
# Create a correlation heatmap
return(pheatmap(mtx,
annotation=HeatmapAnno,
main=paste("Sample Correlation Heatmap:",
Title)))
}Sample QC: Principal Component Analysis
- Checkpoints in PCA: source of variation, sample outlier
QCPCA_fn(TPM, "QC PCA: TPM")
QCPCA_fn(Counts, "QC PCA: Counts")
Sample QC: Sample Correlation Heatmap
- Checkpoints of correlation heatmap: distance between samples, correlation in a group
- Upper: TPM input
- Lower: Counts input
# TPM
QCcorrHeatmap_fn(TPM, "TPM") 
# Counts
QCcorrHeatmap_fn(Counts, "Counts") 
Results & Discussion 2:
- Unsupervised learning-mediated sample QC resulted in highly similar patterns between TPM and Counts
Running DE analysis
# Create a list for TPM and Counts dds
ddsList <- list(TPM=TPM[[1]], Counts=Counts[[1]])
for (x in TvC) {
# Run DESeq()
ddsList[[x]] <- DESeq(ddsList[[x]])
print(resultsNames(ddsList[[x]]))
}## [1] "Intercept" "Group_B_vs_A"
## [1] "Intercept" "Group_B_vs_A"Creating dispersion plots
- Dispersion is important since estimation by DESeq2 algorithm is based on the assumption that genes with similar expression levels have similar dispersion. If an RNA-seq dataset doesn’t satisfy this assumption, use other DE algorithms than DESeq2.
- References: DESeq2 doc "Dispersion plot and fitting alternatives, Harvard Chan Bioinformatics Core workshop
# Plot dispersion
for (x in TvC) {
plotDispEsts(ddsList[[x]],
ylab="Dispersion",
xlab="Mean of Normalized Counts",
main=paste("Dispersion of", x, "Input"))
}
Results & Discussion 3:
- Both inputs gave a goot fit
- No major difference was seen in the dispersion plots between two conditions
Performing shrinkage
- Shrinkage reduces false positives
- Magnitude of shrinkage is affected by dispersion and sample size
- When the argument type is set to “apeglm”, the coef argument is used instead of constrast. In this dataset, you can set coef=Coef where Coef <- resultsNames(ddsList[[1]]).
- When the type is set to “normal”, the argument contrast is set as shown below.
- References: DESeq2 doc “Alternative shrinkage estimators”, Harvard Chan Bioinformatics Core workshop
# Create an empty list for shrunken dds
shr_ddsList <- list(TPM=c(), Counts=c())
for (x in TvC) {
# shrink
shr_ddsList[[x]] <- lfcShrink(ddsList[[x]],
contrast=Contrast, # contrast
type="normal") # is paired with "normal" type
}Extracting log2FoldChange and p-values with or without shrinkage
- The alpha denotes threshold of false discovery rate (FDR) assigned by users.
- In this analysis, the alpha is set to 0.1
# Set FDR threshold
alpha=0.1
# FDR threshold vector
FDRv=c("< 0.1", "> 0.1")
# Set a function cleaning table
Sig_fn <- function(df, Input) {
df <- df %>%
rownames_to_column(var="Gene") %>%
mutate(FDR=ifelse(padj < 0.1 & !is.na(padj),
FDRv[1],
FDRv[2]),
Input=Input)
return(df)
}
# Initialize lists of result tables with (resList) or without (shr_resList) shrinkage
resList <- ddsList
shr_resList <- ddsList
for (x in TvC) {
# Extract results
resList[[x]] <- as.data.frame(results(ddsList[[x]],
contrast=Contrast,
alpha=alpha))
shr_resList[[x]] <- as.data.frame(shr_ddsList[[x]])
# clean the data frame
resList[[x]] <- Sig_fn(resList[[x]], x)
shr_resList[[x]] <- Sig_fn(shr_resList[[x]], x)
}Exploratory data analysis of the extracted log2FoldChange tables
# No shrinkage summary
summary(resList)## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame listhead(resList[['TPM']])## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 11.2413583 0.7644805 0.7163354 1.0672102 2.858769e-01
## 2 A3GALT2 3.7110877 4.2586719 1.6186278 2.6310384 8.512442e-03
## 3 A4GNT 3.8038346 1.2723597 1.1529779 1.1035421 2.697918e-01
## 4 AACSP1 55.7607695 -2.1138938 0.3648669 -5.7936029 6.889223e-09
## 5 AADACL2 0.1641094 0.8714154 4.0804729 0.2135575 8.308922e-01
## 6 AADACL4 0.5008675 -2.5075764 3.7818782 -0.6630506 5.072982e-01
## padj FDR Input
## 1 3.534571e-01 > 0.1 TPM
## 2 1.468366e-02 < 0.1 TPM
## 3 3.361484e-01 > 0.1 TPM
## 4 1.735453e-08 < 0.1 TPM
## 5 NA > 0.1 TPM
## 6 NA > 0.1 TPMhead(resList[['Counts']])## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 10.8066609 0.6941601 0.7188686 0.9656286 3.342301e-01
## 2 A3GALT2 3.8166660 4.2974513 1.6078205 2.6728426 7.521152e-03
## 3 A4GNT 3.8610666 1.2967288 1.1440446 1.1334600 2.570211e-01
## 4 AACSP1 56.5992344 -2.1461615 0.3638336 -5.8987446 3.662778e-09
## 5 AADACL2 0.1700555 0.9200020 4.0804729 0.2254645 8.216179e-01
## 6 AADACL4 0.5019710 -2.4653985 3.7632530 -0.6551243 5.123877e-01
## padj FDR Input
## 1 4.052284e-01 > 0.1 Counts
## 2 1.302526e-02 < 0.1 Counts
## 3 3.212300e-01 > 0.1 Counts
## 4 9.269316e-09 < 0.1 Counts
## 5 NA > 0.1 Counts
## 6 NA > 0.1 Counts# Shrinkage summary
summary(shr_resList)## Length Class Mode
## TPM 9 data.frame list
## Counts 9 data.frame listhead(shr_resList[['TPM']])## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 11.2413583 0.7473037 0.6999647 1.0672102 2.858769e-01
## 2 A3GALT2 3.7110877 3.7983247 1.4378564 2.6310384 8.512442e-03
## 3 A4GNT 3.8038346 1.2010948 1.0856049 1.1035421 2.697918e-01
## 4 AACSP1 55.7607695 -2.1013085 0.3625105 -5.7936029 6.889223e-09
## 5 AADACL2 0.1641094 0.4978588 2.3312198 0.2135575 8.308922e-01
## 6 AADACL4 0.5008675 -1.5560865 2.3037874 -0.6630506 5.072982e-01
## padj FDR Input
## 1 3.534571e-01 > 0.1 TPM
## 2 1.468366e-02 < 0.1 TPM
## 3 3.361484e-01 > 0.1 TPM
## 4 1.735453e-08 < 0.1 TPM
## 5 NA > 0.1 TPM
## 6 NA > 0.1 TPMhead(shr_resList[['Counts']])## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1 A2MP1 10.8066609 0.6783316 0.7023183 0.9656286 3.342301e-01
## 2 A3GALT2 3.8166660 3.8371575 1.4296784 2.6728426 7.521152e-03
## 3 A4GNT 3.8610666 1.2249101 1.0778156 1.1334600 2.570211e-01
## 4 AACSP1 56.5992344 -2.1333964 0.3614677 -5.8987446 3.662778e-09
## 5 AADACL2 0.1700555 0.5246799 2.3270647 0.2254645 8.216179e-01
## 6 AADACL4 0.5019710 -1.5324064 2.2975661 -0.6551243 5.123877e-01
## padj FDR Input
## 1 4.052284e-01 > 0.1 Counts
## 2 1.302526e-02 < 0.1 Counts
## 3 3.212300e-01 > 0.1 Counts
## 4 9.269316e-09 < 0.1 Counts
## 5 NA > 0.1 Counts
## 6 NA > 0.1 CountsExploring mean-difference relationship with MA plots
- Reference: DESeq2 doc “MA-plot”
# Set ylim: has to adjusted by users depending on data
yl <- c(-20, 20)
# Set min log2 fold change of interest
mLog <- c(-1, 1)
# Define a function creating an MA plot
MA_fn <- function(List, Shr) {
MAList <- ddsList
for (i in 1:2) {
MAplot <- ggplot(List[[i]],
aes(x=baseMean,
y=log2FoldChange,
color=FDR)) + geom_point() + scale_x_log10() + theme_bw() + scale_color_manual(values=c("blue", "grey")) + ggtitle(paste("MA plot:", names(List)[i], "Input with", Shr)) + ylim(yl[1], yl[2])+ geom_hline(yintercept=c(mLog[1], mLog[2]), linetype="dashed", color="red")
MAList[[i]] <- MAplot
}
return(MAList)
}
# Create MA plots with or without shrinkage and store in a list
MA <- MA_fn(resList, "No Shrinkage")
shr_MA <- MA_fn(shr_resList, "Shrinkage")Displaying MA plots
- x-axis: expression level (baseMean))
- y-axis: fold change (log2FoldChange)
- Red dashed lines: log2FoldChange = -1 and 1
- Upper: TPM with (right) or without (left) shrinkage
- Lower: Counts with (right) or without (left) shrinkage
# TPM with or without shrinkage
grid.arrange(MA[[1]], shr_MA[[1]], nrow=1)
# TPM with or without shrinkage
grid.arrange(MA[[2]], shr_MA[[2]], nrow=1)
Results & Discussion 4:
- Both inputs gave fairly identical MA plots
- No striking effect was seen by shrinkage
Exploring distribution of false discovery rate (FDR)
- Only the shrunken results are taken for further analysis
- Distribution of adjusted p-val (FDR) was presented
- x-axis: FDR
- y-axis: Number of genes
- Black dashed line: FDR = 0.1
# Combining total data table
res <- rbind(resList[['TPM']], resList[['Counts']])
res$Input <- factor(res$Input, levels=TvC) # TvC=c("TPM", "Counts")
# Create a plot presenting distribution of FDR
ggplot(res,
aes(x=padj, color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ggtitle("Distribution of False Discovery Rate (FDR)") +
xlab("Adjusted P-Value (FDR)") +
ylab("Count") +
geom_vline(xintercept=alpha,
color="black",
linetype="dashed",
size=1) +
scale_x_continuous(breaks=seq(0, 1, by=0.1)) 
Results & Discussion 5:
- Both inputs showed highly similar FDR distribution patterns across the FDR range.
- When the FDR > 0.7, more genes were counted in the TPM input than in the Counts.
Exploring distribution of log2FoldChange by input type
- Black dashed lines: log2FoldChange = -1 and 1
- x-axis: gene expression level (log2FoldChange)
- y-axis: number of genes
ggplot(res[res$FDR == "< 0.1", ], # Subset rows with FDR < alpha
aes(x=log2FoldChange,
color=Input)) +
geom_density(size=1, aes(y=..count..)) +
theme_bw() +
ylab("Count") +
geom_vline(xintercept=c(mLog[1], mLog[2]),
color="black",
linetype="dashed", size=1) +
ggtitle("Distribution of Log2 Folds by Input Type") +
xlim(xlim[1], xlim[2])
Results & Discussion 7:
- Both inputs showed nearly identical fold change distribution
Exploring expression profiling with normalized count data
- Normalized count matrices are extracted from dds objects and filtered with thresholds set at FDR and log2FoldChange
- The heatmaps display z-scores of the normalized counts
- lowfdrList: a list of matrices filtered by FDR < alpha
- highfoldList: a list of matrices filtered by FDR < alpha AND absolute log2FoldChange > user’s minimum threshold (mLog)
- In this analysis, mLog = 1
- References: Harvard Chan Bioinformatics Core workshop, DESeq2 doc “Heatmap of the count matrix”
# Initialize a list
lowfdrList <- ddsList # A list for normalized counts matrix with FDR below alpha
highfoldList <- ddsList # A list for normalized counts with log2foldchange over minLog
for (x in TvC) {
# Create filtering vectors with alpha and log2foldchange
BelowAlpha <- resList[[x]]$FDR == FDRv[1]
overmLog <- abs(resList[[x]]$log2FoldChange) > mLog[2] # mLog has been set to c(-1, 1) previously
# Extract transformed counts from vsd objects (TPM[['vsd']] and Counts[['vsd']])
if (x == "TPM") {
normCounts <- counts(TPM[['dds']], normalized=T)
} else {
normCounts <- counts(Counts[['dds']], normalized=T)
}
# Update the normalized count matrix with FDR below alpha
lowfdrList[[x]] <- normCounts[BelowAlpha, ]
highfoldList[[x]] <- normCounts[BelowAlpha & overmLog, ]
summary(lowfdrList[[x]])
summary(highfoldList[[x]])
}
# Initialize map lists
lowfdrMap <- ddsList
highfoldMap <- ddsList
# Set a function creating a heatmap
ProfileHeatmap_fn <- function(inputmatrix, Title1, Title2, Title3=NULL) {
as.ggplot(pheatmap(inputmatrix,
annotation=HeatmapAnno,
scale="row", # presents z-score instead of counts
show_rownames=F,
main=paste("Transcription Profiles with",
Title1,
"input and",
Title2,
alpha,
Title3)))
}
# Create and save heatmaps
for (x in TvC) {
lowfdrMap[[x]] <- ProfileHeatmap_fn(lowfdrList[[x]],
Title1=x,
Title2="FDR <")
highfoldMap[[x]] <- ProfileHeatmap_fn(highfoldList[[x]],
Title1=x,
Title2="FDR <",
Title3=paste("+ Absolte Log2 Fold Change >", mLog[2]))
}
Results & Discussion 8:
- No major difference is seen bewteen the TPM and Counts
NA statistics: zero count genes & outlier genes
When NAs appear in
- log2FoldChange: zero counts in all samples
- padj: too little information
- pval & padj: at least one replicate was an outlier
# Count number of NA genes
type=c("Zero Counts", "Outliers", "Total NA Genes")
NAstat <- res %>%
group_by(Input) %>%
summarize(zero=sum(is.na(log2FoldChange)),
outlier=sum(is.na(pvalue) & is.na(padj))) %>%
mutate(total=sum(zero, outlier)) %>%
gather(Type, Number, -Input) %>%
mutate(Type=factor(case_when(Type == "zero" ~ type[1],
Type == "outlier" ~ type[2],
Type == "total" ~ type[3]),
levels=type))
# Plot number of NA genes
ggplot(NAstat, aes(x=Type, y=Number, group=Input, fill=Input, label=Number)) +
geom_bar(stat="identity", position="dodge") +
theme_bw() +
geom_text(position=position_dodge(width=1), vjust=1.5) +
ggtitle("Number of NA Genes") +
ylab("Number of Genes")
Results & Discussion 9:
- Both inputs produced no major difference in the number of NA genes
Ranking DEGs with the TPM and original count inputs
FDRrankList: ranked by FDR
lfcList: ranked by absolute fold change
UPlfcList: ranked by magnitude of fold change increase
DOWNlfcList: ranked by manitude of fold change decrease
# Create a new list having DE table with FDR below alpha
lowfdr_resList <- resList
for (x in TvC) {
lowfdr_resList[[x]] <- filter(resList[[x]],
FDR == FDRv[1]) %>%
as.data.table()
}
# Initialize new lists in order to store rank-updated result DE tables
FDRrankList <- lowfdr_resList
lfcList <- lowfdr_resList
UPlfcList <- lowfdr_resList
DOWNlfcList <- lowfdr_resList
# Set a function creating a column for the rank
Ranking_fn <- function(x) {mutate(x, Rank=1:nrow(x))}
for (x in TvC) {
# Rearrange genes with FDR
FDRrankList[[x]] <- lowfdr_resList[[x]][order(padj),] %>%
Ranking_fn()
# Rearrange genew with absolute log2FoldChange
lfcList[[x]] <- lowfdr_resList[[x]][order(-abs(log2FoldChange)),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (decreasing order)
UPlfcList[[x]] <- lowfdr_resList[[x]][order(-log2FoldChange),] %>%
Ranking_fn()
# Rearrange genes with log2FoldChange (increasing order)
DOWNlfcList[[x]] <- lowfdr_resList[[x]][order(log2FoldChange),] %>%
Ranking_fn()
}
# Explore the ranks
print(c(FDRrankList, lfcList, UPlfcList, DOWNlfcList))## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: ACBD7 3168.682742 2.9056054 0.06424347 45.228029 0.00000000
## 2: AFP 32584.595366 -11.0421341 0.15395409 -71.723551 0.00000000
## 3: AKAP11 13905.676097 2.4997058 0.03666677 68.173600 0.00000000
## 4: AMIGO1 3568.001961 4.4114268 0.06849514 64.404960 0.00000000
## 5: AMY2B 1557.113021 3.6430792 0.08662900 42.053805 0.00000000
## ---
## 4726: SALL4P1 1.454384 -4.0531916 2.21570166 -1.829304 0.06735410
## 4727: GDF9 70.162940 0.5210032 0.28544731 1.825217 0.06796832
## 4728: RPL12P10 3.588481 -2.6723591 1.46472861 -1.824474 0.06808046
## 4729: RPL34P18 9.116960 -1.4682061 0.80491820 -1.824044 0.06814545
## 4730: IMPDH1P10 12.307808 -7.1302067 3.91051047 -1.823344 0.06825127
## padj FDR Input Rank
## 1: 0.00000000 < 0.1 TPM 1
## 2: 0.00000000 < 0.1 TPM 2
## 3: 0.00000000 < 0.1 TPM 3
## 4: 0.00000000 < 0.1 TPM 4
## 5: 0.00000000 < 0.1 TPM 5
## ---
## 4726: 0.09876745 < 0.1 TPM 4726
## 4727: 0.09964469 < 0.1 TPM 4727
## 4728: 0.09978798 < 0.1 TPM 4728
## 4729: 0.09986212 < 0.1 TPM 4729
## 4730: 0.09999605 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: ABCD2 1069.176400 5.560740 0.14812494 37.540872 0.00000000
## 2: ACBD7 3223.871084 2.902042 0.06325036 45.881828 0.00000000
## 3: AFP 32279.951238 -11.043880 0.14081349 -78.429135 0.00000000
## 4: AKAP11 14142.443732 2.494980 0.03650368 68.348719 0.00000000
## 5: AMIGO1 3630.556248 4.406517 0.06828149 64.534572 0.00000000
## ---
## 4730: SUGT1P2 10.591190 -1.285347 0.70388149 -1.826085 0.06783752
## 4731: MTND2P5 4.330142 2.185290 1.19716309 1.825390 0.06794219
## 4732: RPSAP21 2.057284 -3.437381 1.88421063 -1.824308 0.06810546
## 4733: OR52W1 3.905048 2.282493 1.25138546 1.823973 0.06815620
## 4734: SALL4P1 1.468341 -4.021370 2.20511483 -1.823656 0.06820416
## padj FDR Input Rank
## 1: 0.00000000 < 0.1 Counts 1
## 2: 0.00000000 < 0.1 Counts 2
## 3: 0.00000000 < 0.1 Counts 3
## 4: 0.00000000 < 0.1 Counts 4
## 5: 0.00000000 < 0.1 Counts 5
## ---
## 4730: 0.09937551 < 0.1 Counts 4730
## 4731: 0.09950780 < 0.1 Counts 4731
## 4732: 0.09972586 < 0.1 Counts 4732
## 4733: 0.09977907 < 0.1 Counts 4733
## 4734: 0.09982819 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1935.1129 -14.42718266 1.18067729 -12.219412 2.448482e-34
## 2: GAGE12C 665.4894 -12.88787310 1.18254611 -10.898411 1.172881e-27
## 3: GAGE13 612.2774 -12.76494282 1.18202949 -10.799175 3.473116e-27
## 4: GAGE12F 2805.4936 -12.51900730 0.83651932 -14.965593 1.232215e-50
## 5: UTP14C 531.7746 12.47219039 1.18872293 10.492092 9.392533e-26
## ---
## 4726: GM2A 5597.9316 0.10758757 0.04117832 2.612724 8.982382e-03
## 4727: CD2BP2 6723.4716 -0.10030164 0.04900183 -2.046896 4.066833e-02
## 4728: UQCRFS1 5731.6147 -0.09801823 0.04466793 -2.194376 2.820840e-02
## 4729: PRPF6 8832.8778 -0.09437979 0.04835202 -1.951931 5.094644e-02
## 4730: SELENOP 20675.0072 -0.08968953 0.04428648 -2.025212 4.284557e-02
## padj FDR Input Rank
## 1: 1.092594e-33 < 0.1 TPM 1
## 2: 4.687467e-27 < 0.1 TPM 2
## 3: 1.378505e-26 < 0.1 TPM 3
## 4: 6.681727e-50 < 0.1 TPM 4
## 5: 3.642432e-25 < 0.1 TPM 5
## ---
## 4726: 1.543081e-02 < 0.1 TPM 4726
## 4727: 6.294319e-02 < 0.1 TPM 4727
## 4728: 4.503207e-02 < 0.1 TPM 4728
## 4729: 7.715446e-02 < 0.1 TPM 4729
## 4730: 6.601152e-02 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1975.0482 -14.41206728 1.18066330 -12.206755 2.860740e-34
## 2: GAGE12C 679.3481 -12.87323165 1.18248182 -10.886621 1.335017e-27
## 3: GAGE13 624.8986 -12.74994147 1.18192899 -10.787401 3.947987e-27
## 4: GAGE12F 2863.4359 -12.51268995 0.83641254 -14.959950 1.341265e-50
## 5: UTP14C 541.2130 12.50175251 1.18864170 10.517680 7.161562e-26
## ---
## 4730: CD2BP2 6800.0952 -0.10470103 0.04950383 -2.115009 3.442919e-02
## 4731: GM2A 5697.3428 0.10317346 0.04064652 2.538310 1.113893e-02
## 4732: UQCRFS1 5839.0452 -0.10252336 0.04427907 -2.315391 2.059153e-02
## 4733: PRPF6 8994.9911 -0.09892368 0.04794184 -2.063410 3.907370e-02
## 4734: SELENOP 21049.9178 -0.09447398 0.04403219 -2.145566 3.190761e-02
## padj FDR Input Rank
## 1: 1.266586e-33 < 0.1 Counts 1
## 2: 5.304090e-27 < 0.1 Counts 2
## 3: 1.554295e-26 < 0.1 Counts 3
## 4: 7.238026e-50 < 0.1 Counts 4
## 5: 2.758336e-25 < 0.1 Counts 5
## ---
## 4730: 5.388748e-02 < 0.1 Counts 4730
## 4731: 1.887083e-02 < 0.1 Counts 4731
## 4732: 3.361892e-02 < 0.1 Counts 4732
## 4733: 6.051445e-02 < 0.1 Counts 4733
## 4734: 5.024724e-02 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: UTP14C 531.7746 12.47219 1.1887229 10.49209 9.392533e-26
## 2: NEUROD6 914.0304 12.29159 1.1806402 10.41096 2.209892e-25
## 3: CARTPT 444.8165 12.21474 1.1829535 10.32563 5.396390e-25
## 4: HBB 2497.0069 11.87435 0.7237799 16.40603 1.731676e-60
## 5: SLC17A6 641.0675 11.78002 1.1827184 9.96012 2.277885e-23
## ---
## 4726: BPIFB2 1748.3190 -12.42289 1.0233841 -12.13902 6.560335e-34
## 4727: GAGE12F 2805.4936 -12.51901 0.8365193 -14.96559 1.232215e-50
## 4728: GAGE13 612.2774 -12.76494 1.1820295 -10.79917 3.473116e-27
## 4729: GAGE12C 665.4894 -12.88787 1.1825461 -10.89841 1.172881e-27
## 4730: RHOXF2B 1935.1129 -14.42718 1.1806773 -12.21941 2.448482e-34
## padj FDR Input Rank
## 1: 3.642432e-25 < 0.1 TPM 1
## 2: 8.512813e-25 < 0.1 TPM 2
## 3: 2.060440e-24 < 0.1 TPM 3
## 4: 1.031859e-59 < 0.1 TPM 4
## 5: 8.459668e-23 < 0.1 TPM 5
## ---
## 4726: 2.904992e-33 < 0.1 TPM 4726
## 4727: 6.681727e-50 < 0.1 TPM 4727
## 4728: 1.378505e-26 < 0.1 TPM 4728
## 4729: 4.687467e-27 < 0.1 TPM 4729
## 4730: 1.092594e-33 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: UTP14C 541.2130 12.50175 1.1886417 10.517680 7.161562e-26
## 2: NEUROD6 929.8364 12.30527 1.1805725 10.423137 1.944330e-25
## 3: CARTPT 452.4859 12.24355 1.1828268 10.351089 4.137522e-25
## 4: HBB 2400.0112 12.13599 0.8113413 14.957939 1.382416e-50
## 5: SLC17A6 652.2483 11.81567 1.1825758 9.991472 1.660961e-23
## ---
## 4730: BPIFB2 1782.0106 -12.45487 1.0233020 -12.171251 4.422530e-34
## 4731: GAGE12F 2863.4359 -12.51269 0.8364125 -14.959950 1.341265e-50
## 4732: GAGE13 624.8986 -12.74994 1.1819290 -10.787401 3.947987e-27
## 4733: GAGE12C 679.3481 -12.87323 1.1824818 -10.886621 1.335017e-27
## 4734: RHOXF2B 1975.0482 -14.41207 1.1806633 -12.206755 2.860740e-34
## padj FDR Input Rank
## 1: 2.758336e-25 < 0.1 Counts 1
## 2: 7.430922e-25 < 0.1 Counts 2
## 3: 1.574349e-24 < 0.1 Counts 3
## 4: 7.454287e-50 < 0.1 Counts 4
## 5: 6.131486e-23 < 0.1 Counts 5
## ---
## 4730: 1.953073e-33 < 0.1 Counts 4730
## 4731: 7.238026e-50 < 0.1 Counts 4731
## 4732: 1.554295e-26 < 0.1 Counts 4732
## 4733: 5.304090e-27 < 0.1 Counts 4733
## 4734: 1.266586e-33 < 0.1 Counts 4734
##
## $TPM
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1935.1129 -14.42718 1.1806773 -12.21941 2.448482e-34
## 2: GAGE12C 665.4894 -12.88787 1.1825461 -10.89841 1.172881e-27
## 3: GAGE13 612.2774 -12.76494 1.1820295 -10.79917 3.473116e-27
## 4: GAGE12F 2805.4936 -12.51901 0.8365193 -14.96559 1.232215e-50
## 5: BPIFB2 1748.3190 -12.42289 1.0233841 -12.13902 6.560335e-34
## ---
## 4726: SLC17A6 641.0675 11.78002 1.1827184 9.96012 2.277885e-23
## 4727: HBB 2497.0069 11.87435 0.7237799 16.40603 1.731676e-60
## 4728: CARTPT 444.8165 12.21474 1.1829535 10.32563 5.396390e-25
## 4729: NEUROD6 914.0304 12.29159 1.1806402 10.41096 2.209892e-25
## 4730: UTP14C 531.7746 12.47219 1.1887229 10.49209 9.392533e-26
## padj FDR Input Rank
## 1: 1.092594e-33 < 0.1 TPM 1
## 2: 4.687467e-27 < 0.1 TPM 2
## 3: 1.378505e-26 < 0.1 TPM 3
## 4: 6.681727e-50 < 0.1 TPM 4
## 5: 2.904992e-33 < 0.1 TPM 5
## ---
## 4726: 8.459668e-23 < 0.1 TPM 4726
## 4727: 1.031859e-59 < 0.1 TPM 4727
## 4728: 2.060440e-24 < 0.1 TPM 4728
## 4729: 8.512813e-25 < 0.1 TPM 4729
## 4730: 3.642432e-25 < 0.1 TPM 4730
##
## $Counts
## Gene baseMean log2FoldChange lfcSE stat pvalue
## 1: RHOXF2B 1975.0482 -14.41207 1.1806633 -12.206755 2.860740e-34
## 2: GAGE12C 679.3481 -12.87323 1.1824818 -10.886621 1.335017e-27
## 3: GAGE13 624.8986 -12.74994 1.1819290 -10.787401 3.947987e-27
## 4: GAGE12F 2863.4359 -12.51269 0.8364125 -14.959950 1.341265e-50
## 5: BPIFB2 1782.0106 -12.45487 1.0233020 -12.171251 4.422530e-34
## ---
## 4730: SLC17A6 652.2483 11.81567 1.1825758 9.991472 1.660961e-23
## 4731: HBB 2400.0112 12.13599 0.8113413 14.957939 1.382416e-50
## 4732: CARTPT 452.4859 12.24355 1.1828268 10.351089 4.137522e-25
## 4733: NEUROD6 929.8364 12.30527 1.1805725 10.423137 1.944330e-25
## 4734: UTP14C 541.2130 12.50175 1.1886417 10.517680 7.161562e-26
## padj FDR Input Rank
## 1: 1.266586e-33 < 0.1 Counts 1
## 2: 5.304090e-27 < 0.1 Counts 2
## 3: 1.554295e-26 < 0.1 Counts 3
## 4: 7.238026e-50 < 0.1 Counts 4
## 5: 1.953073e-33 < 0.1 Counts 5
## ---
## 4730: 6.131486e-23 < 0.1 Counts 4730
## 4731: 7.454287e-50 < 0.1 Counts 4731
## 4732: 1.574349e-24 < 0.1 Counts 4732
## 4733: 7.430922e-25 < 0.1 Counts 4733
## 4734: 2.758336e-25 < 0.1 Counts 4734Comparing DEG ranks between TPM- and Counts-inputted DE analysis
# Set a function rebuilding DE tables with gene ranks
combineTable_fn <- function(List){
# Select columns and join the data frames by gene
full_join(List[[TvC[1]]][,.(Gene, Input, Rank, baseMean)],
List[[TvC[2]]][,.(Gene, Input, Rank, baseMean)], by="Gene") %>%
# Add columns assining gene expression levels, rank differences, and mean ranks
mutate(logMeanExpression=log(baseMean.x+baseMean.y/2),
RankDiff=Rank.x-Rank.y,
MeanRank=(Rank.x+Rank.y)/2)
}
# Explore outputs of the function
head(combineTable_fn(FDRrankList))## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y logMeanExpression
## 1: ACBD7 TPM 1 3168.683 Counts 2 3223.871 8.472325
## 2: AFP TPM 2 32584.595 Counts 3 32279.951 10.793939
## 3: AKAP11 TPM 3 13905.676 Counts 4 14142.444 9.951177
## 4: AMIGO1 TPM 4 3568.002 Counts 5 3630.556 8.591053
## 5: AMY2B TPM 5 1557.113 Counts 6 1581.779 7.761320
## 6: ANP32B TPM 6 10924.599 Counts 7 11143.790 9.710903
## RankDiff MeanRank
## 1: -1 1.5
## 2: -1 2.5
## 3: -1 3.5
## 4: -1 4.5
## 5: -1 5.5
## 6: -1 6.5dim(combineTable_fn(FDRrankList))## [1] 4761 10tail(combineTable_fn(FDRrankList))## Gene Input.x Rank.x baseMean.x Input.y Rank.y baseMean.y
## 1: INSL6 <NA> NA NA Counts 4723 8.692921
## 2: ASS1P10 <NA> NA NA Counts 4726 9.072376
## 3: RPS3AP30 <NA> NA NA Counts 4729 1.481576
## 4: SUGT1P2 <NA> NA NA Counts 4730 10.591190
## 5: MTND2P5 <NA> NA NA Counts 4731 4.330142
## 6: OR52W1 <NA> NA NA Counts 4733 3.905048
## logMeanExpression RankDiff MeanRank
## 1: NA NA NA
## 2: NA NA NA
## 3: NA NA NA
## 4: NA NA NA
## 5: NA NA NA
## 6: NA NA NA# Set a function plotting gene ranks between TPM- (x-axis) and Counts-Inputs (y-axis)
compareRanks_fn <- function(df, rankedby) {
ggplot(df,
aes(x=Rank.x,
y=Rank.y,
color=logMeanExpression)) +
geom_point(alpha=0.5) +
theme_bw() +
xlab("Rank with TPM") +
ylab("Rank with Counts") +
geom_abline(slope=1, color="black", size=0.5) +
ggtitle(paste(rankedby, "Ranking with TPM vs Counts Inputs")) +
scale_color_gradient(low="blue", high="red")
}
# Set a function plotting the rank difference over the gene expression level
RankdiffOverMean_fn <- function(df, rankedby) {
ggplot(df,
aes(x=logMeanExpression,
y=RankDiff,
color=MeanRank)) +
geom_point(alpha=0.5) +
theme_bw() +
ylab("Rank Difference (TPM - Counts)") +
ggtitle(paste("Rank Difference Inputs (TPM - Counts) in", rankedby)) +
geom_hline(yintercept=0, color="black", size=0.5) +
scale_color_gradient(low="blue", high="red")
}Visualizing DEG ranks I: TPM- vs Counts-input
- x-axis: rank with TPM input
- y-axis: rank with Counts input
- Black diagonal lines: rank with TPM = rank with Counts
- Dot color: gene expression level (log-baseMean)
# Ranked by FDR
compareRanks_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
compareRanks_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
compareRanks_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
compareRanks_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
Results & Discussion 10:
- FDR ranking (top): Genes with high expression levels (red dots) were ranked similarly between two input types which is close to the black diagonal line. Ranking dissimilarity increased in lower ranking ranges (250-500).
- Absolute fold ranking (2nd top): Genes were ranked higher as their expression levels were lower (blue to purple dots). Genes with high expression levels (red dots) showed less difference in the rankings which is close to the black diagonal line.
- Decrease/increase fold ranking (bottom/2nd botton): Highly expressed genes (red dots) showed middle-rankings (150-350). Low-ranking genes tended to be ranked higher with TPM than with Counts as the slope was higher than 1 which is indicated by the black diagonal line.
Visualizing DEG ranks II: Relationship between gene expression level and rank difference
- x-axis: expression level (log-baseMean)
- y-axis: rank difference (rank with TPM - rank with Counts)
- Black horizontal lines: rank with TPM = rank with Counts
- Dot color: average rank
# Ranked by FDR
RankdiffOverMean_fn(combineTable_fn(FDRrankList),
"FDR")
# Ranked by absolute fold change
RankdiffOverMean_fn(combineTable_fn(lfcList),
"Absolute Log2FoldChange")
# Ranked by magnitude of positive fold change
RankdiffOverMean_fn(combineTable_fn(UPlfcList),
"Log2FoldChange (Increased)")
# Ranked by magnitude of negative fold change
RankdiffOverMean_fn(combineTable_fn(DOWNlfcList),
"Log2FoldChange (Decreased)")
Results & Discussion 11:
- Difference in FDR ranking (top): High-ranking genes (blue dots) showed less difference. Genes with high expression levels (higher x-values) showed less difference.
- Difference in absolute/increase fold ranking (2nd top/2nd bottom): Low mean ranks (red dots) seemed to be coupled with higher ranking with TPM (negative y-values) than that with Counts while high mean ranks (purple to blue dots) seemed vice versa or random.
- Difference in decrease fold ranking (bottom): Low-to-middle ranking genes (red-purple dots) tended to have higher ranking with TPM (negative y-values) whereas high ranking genes (blue dots) showed higher ranking with Counts (positive y-values)
Summarizing up/down DEGs with an upset plot
- Calculate the number of genes
# Clean data to generate an upset plot
res <- res %>%
# Filter genes with valid padj
filter(!is.na(padj)) %>%
# Detect genes which are up/down/unchanged change patterns in either TPM and Counts inputs
mutate(Up=ifelse(FDR == FDRv[1] & log2FoldChange > 0, Gene, ""), # What are upregulated genes?
Down=ifelse(FDR == FDRv[1] & log2FoldChange < 0, Gene, ""), # What are downregulated genes?
Unchanged=ifelse(FDR == FDRv[2], Gene, ""), # What are unchanged genes?
TPM_Input=ifelse(Input == "TPM", Gene, ""), # What are the genes from TPM input?
Counts_Input=ifelse(Input == "Counts", Gene, "")) # What are the genes from Counts input?
# Create a list storing groups of interest
upsetInput <- list(Up=res$Up,
Down=res$Down,
Unchanged=res$Unchanged,
TPM_Input=res$TPM,
Counts_Input=res$Counts)
# Create the upset plot
upset(fromList(upsetInput),
sets.x.label="Number of Genes per Group",
order.by="freq") 
Session Info
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/mira/miniconda3/envs/r/lib/libopenblasp-r0.3.10.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] AnnotationDbi_1.52.0 UpSetR_1.4.0
## [3] ggrepel_0.8.2 ggplotify_0.0.5
## [5] gridExtra_2.3 pheatmap_1.0.12
## [7] DESeq2_1.30.0 SummarizedExperiment_1.20.0
## [9] Biobase_2.50.0 MatrixGenerics_1.2.0
## [11] matrixStats_0.57.0 GenomicRanges_1.42.0
## [13] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [15] S4Vectors_0.28.0 tximport_1.18.0
## [17] forcats_0.5.0 stringr_1.4.0
## [19] dplyr_1.0.2 purrr_0.3.4
## [21] readr_1.4.0 tidyr_1.1.2
## [23] tibble_3.0.4 ggplot2_3.3.2
## [25] tidyverse_1.3.0 AnnotationHub_2.22.0
## [27] BiocFileCache_1.14.0 dbplyr_2.0.0
## [29] BiocGenerics_0.36.0 rmarkdown_2.5
## [31] data.table_1.13.2
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1
## [3] XVector_0.30.0 fs_1.5.0
## [5] rstudioapi_0.13 farver_2.0.3
## [7] bit64_4.0.5 interactiveDisplayBase_1.28.0
## [9] fansi_0.4.1 lubridate_1.7.9.2
## [11] xml2_1.3.2 splines_4.0.3
## [13] geneplotter_1.68.0 knitr_1.30
## [15] jsonlite_1.7.1 broom_0.7.2
## [17] annotate_1.68.0 shiny_1.5.0
## [19] BiocManager_1.30.10 compiler_4.0.3
## [21] httr_1.4.2 rvcheck_0.1.8
## [23] backports_1.2.0 assertthat_0.2.1
## [25] Matrix_1.2-18 fastmap_1.0.1
## [27] cli_2.2.0 later_1.1.0.1
## [29] htmltools_0.5.0 tools_4.0.3
## [31] gtable_0.3.0 glue_1.4.2
## [33] GenomeInfoDbData_1.2.4 rappdirs_0.3.1
## [35] Rcpp_1.0.5 cellranger_1.1.0
## [37] vctrs_0.3.5 xfun_0.19
## [39] ps_1.5.0 rvest_0.3.6
## [41] mime_0.9 lifecycle_0.2.0
## [43] XML_3.99-0.3 zlibbioc_1.36.0
## [45] scales_1.1.1 hms_0.5.3
## [47] promises_1.1.1 RColorBrewer_1.1-2
## [49] yaml_2.2.1 curl_4.3
## [51] memoise_1.1.0 stringi_1.5.3
## [53] RSQLite_2.2.1 BiocVersion_3.12.0
## [55] genefilter_1.72.0 BiocParallel_1.24.0
## [57] rlang_0.4.9 pkgconfig_2.0.3
## [59] bitops_1.0-6 evaluate_0.14
## [61] lattice_0.20-41 labeling_0.4.2
## [63] bit_4.0.4 tidyselect_1.1.0
## [65] plyr_1.8.6 magrittr_2.0.1
## [67] R6_2.5.0 generics_0.1.0
## [69] DelayedArray_0.16.0 DBI_1.1.0
## [71] pillar_1.4.7 haven_2.3.1
## [73] withr_2.3.0 survival_3.2-7
## [75] RCurl_1.98-1.2 modelr_0.1.8
## [77] crayon_1.3.4 utf8_1.1.4
## [79] locfit_1.5-9.4 grid_4.0.3
## [81] readxl_1.3.1 blob_1.2.1
## [83] reprex_0.3.0 digest_0.6.27
## [85] xtable_1.8-4 httpuv_1.5.4
## [87] gridGraphics_0.5-0 munsell_0.5.0